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Abstract 


By means of a simple structural problem, we bring into focus an important 
requirement often overlooked in practice on the basis functions used in Rayleigh - 
Ritz-Galerkin type methods. The problem of the static deformation of a uniformly 
loaded beam is solved variationally by expanding the beam displacement in a 
Fourier Cosine series. The potential energy functional is rendered stationary 
subject to the geometric boundary conditions. It is demonstrated that the 
variational approach does not converge to the true solution. The object of the 
paper Is to resolve this paradox, and in so doing, indicate the practical 
implications of norms and completeness in an appropriate Inner product space. 

Introduction 


Virtually all flexible multibody codes in use today are based upon some 
variational principle of mechanics. The most common of these being Hamilton's 
Principle and its discretized version-Lagrange' s Equations. Regardless of the 
particular label attached to the technique (e.g. , assumed modes method, Ritz- 
Galerkin) , the problem reduces to rendering stationary a certain definite 
integral with respect to a sufficiently regular family of functions subject to 
certain "geometric" boundary conditions. In practice the basis functions are 
generated using a general purpose Finite Element program and may be subject to 
further manipulation such as modal synthesis before being incorporated Into the 
multibody program. It is assumed by many analysts that the basis functions used 
are to a certain extent arbitrary. It is argued that so long as they are members 
of an infinite family of orthogonal functions and satisfy the geometric boundary 
conditions, convergence of the dynamic response is guaranteed. 

It is well-known in the mathematical theory of variational methods that the 
basis functions must be complete in an appropriate inner product space. Our 
problem demonstrates the importance of this somewhat subtle criterion. It is 
extremely significant that application of the variational principle using an 
inappropriate set of basis functions is convergent, but to an erroneous result. 
This would be difficult if not impossible to identify in a typical spacecraft 
application 

Problem Statement and Exact Solution 

A cantilever beam constant bending stiffness El and unit length is acted on by 
a uniform distributed loading. If the load per unit length is q , and we write p 
for the ratio q/EI , then the static deflection W(x) is governed by 
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- P, W(0) - W'(0) - W"(l) - W'"(l) - 0 
dx 


With constant p the solution is a polynomial of degree four: 



For future reference we develop this solution in a Fourier cosine series, 
extending it as an even function to (- 1,0): 
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Variational Solution 


A variational principle equivalent to the boundary value problem (1) is the 
principle of minimum total potential energy: Find the function W(x) satisfying 

W(0) — W' (0) * 0 (and sufficiently regular) that minimizes the total potential 
energy El P(W), where 


P(W) - 
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We apply the Rayleigh-Ritz method to the energy functional in (4) . The set of 
basis functions will be {1, cos ?rx, cos 2 jtx,...}, which is complete over the 
interval 0<x<i. The variational solution is then a cosine expansion 
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which automatically satisfies the geometric boundary condition W' (0) - 0 (as well 
as the natural boundary condition W' * 9 (1) - 0) . Substituting the assumed 
expansion (5) into the integral (4) , the orthogonality of the cosines gives the 
value J for the total potential energy: 
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We calculate the coefficients a 0 , a lf ... by rendering J stationary subject to the 
remaining geometric boundary condition: W(0) — 0, or Za„ — 0. 


Following the standard method we introduce a Lagrange multiplier A and build the 
constraint into the auxiliary function 
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The equation 3L/3a 0 - 0 yields A - p. Then dh/da^ - 0 gives 
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The geometric constraint then yields 
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When the coefficients (7) , (8) are inserted into (5) we obtain the variational 
solution 
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This can also be summed exactly, and it is again a fourth degree polynomial - but 
a different one! 


W*(x) 




( 10 ) 


Comparing with (2), the difference W-W* is px 2 /12. Only geometric boundary 
conditions were enforced in computing W*. One expects that the natural boundary 

conditions will automatically be satisfied. Note however that (W*)" p/6 at 

x - 1, and not zero. A plot of W (exact) and W* (Ritz) is given in Figure 1. 
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The paper started by expanding the exact solution W into the cosine series 
(3). That series has coefficients of order 1/n 2 , and it converges to W. But 
computing the bending energy in these two forms of the solution, polynomial and 
Fourier, leads to a disturbing result: 

For the polynomial, 



is finite . 


For the series, 



2 2 
(W M ) dx *■ p 



-2 

4 4 
n * 


+ 


(-l) n ~ 

, 2 2 
3n tt J 


is infinite . 


The orthogonality of the cosines produces this sum of positive terms, roughly 
n^ 2 , and the sum does not converge. 
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Resolution of the Paradox 


When the geometric boundary conditions are linear and homogeneous, it is 
standard practice to choose the basis functions such that each individually 
satisfies them. Since our basis functions of cosines do not all vanish at x=0 , 
one might be tempted to point to this fact as the root of the problem. However, 
our enforcement of the geometric boundary condition through a Lagrange multiplier 
is perfectly legitimate. There is nothing in the theory that forces us to 
satisfy the geometric boundary conditions by each basis function; only the 
resultant linear combination must satisfy them. (See the section "Numeirical 
Results" for substantiation of this statement.) The key to resolving the paradox 
lies in a loose statement preceding eq. (5), regarding the "completeness" of the 
cosines. The word "complete" in itself has no meaning. We have to identify the 
space of admissible functions, as well as its norm (the measure of distance in 
that space) , before it can be claimed that a set of trial functions is complete - 
in other words, before we can say that the combinations of the trial functions 
can approximate with arbitrary accuracy any admissible function. To discuss 
accuracy we need a norm. 

The most common measure of distance is the L z norm — the square root of 
Jf z dx. The function f need not be continuous; step functions present absolutely 
no difficulty. The space contains functions much worse than that, although a 
delta function has infinite length and is not allowed. In the L z norm the set 
(1, cos irx.cos is complete (on the interval 0<x<l) and this is a 

cornerstone of Fourier analysis. The cosines are also complete in the L p spaces, 
with norm (/If^dx) ' p and p>l , but L z is special: it is associated with an 
inner product. That makes it a Hilbert space, in which (f, g) - Jfg dx matches 
the norm: (f,f) agrees with ||f|| 2 . The L 2 space admits angles, and 

orthogonality, and all the geometry of ordinary Euclidean space. But it is not 
the only Hilbert space, nor is it necessarily the right one. 

Our fourth-order problem comes with its own norm and inner product and space 
of admissible functions. We look there for the resolution of the paradox; we 
have to work with the right space. The norm comes directly from the bending 
energy : 

1 

||w|| 2 - j (W"(x)) 2 dx 
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The inner product is determined by the norm: 


(W,V) 



W" (x)V" (x)dx 


The admissible functions are also determined: Their norm must be finite and they 

must satisfy the essential boundary conditions. Thus W is admissible if 

1 

(W") 2 dx < « and W(0)=0 and W'(0)— 0 
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Such a function W comes from integrating twice a function in L 2 : 


W(x) 



f(t)dt dy 


Notice that W is not required to satisfy the natural conditions W" (1)=0 and 
W'"(i)-0. We could not make that requirement and still have a complete space. 
Functions that satisfy these extra conditions can come arbitrarily close 
(measured by the norm) to functions that don't. The process of completion wipes 
out the natural conditions as a requirement on admissible W. In Reference 1 the 
second author referred to the admissible space as H E Z - the Hilbert space of 
functions that have two derivatives (in L 2 - H°) and that satisfy the essential 
boundary conditions. Remark 3 will justify more fully the choice of bending 
energy - the second- degree term in the total potential energy - as the norm. 


That finishes the functional analysis. It was needed in order to ask the 
right question: Are the cosines complete in the space of admissible functions? 

We suspect that the answer must be no. 

Apart from boundary conditions, we are asking whether combinations of the 
cosines (including cos 0-1) can come arbitrarily close to W. We know they can do 
so in the ordinary L 2 norm, and the Fourier expansion (3) does it explicitly. 

The question is whether the cosines can come arbitrarily close in the 'second- * 
derivative norm 1 . Equivalently, the second derivatives of our set of cosines 
must come close, in the ordinary sense to W tt . But the second derivatives are 


{0, -7T 2 cos ttx , -4?r 2 cos 2 jtx,...} 


You see the problem. We are missing the constant term! We cannot approximate 
W n -2 with the functions we have left. In other words, we cannot come close to 
W-x 2 with our original set of trial functions. The cosines were not complete, 
but if this additional trial function x 2 is included, the set is complete. That 
would add the constant f unc tion to the list of second derivatives. So it was no 
accident that the discrepancy between W and W* was a multiple of x 2 . 


Before drawing a final conclusion about W* we add four observations. 

1. The new function x 2 not only completes the set, it is orthogonal to the 
original cosines. The Inner product is 


/ 


cos nirx ■ 2 dx - 0 
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2. The original cosines were not in H E 2 - they were not really admissible - 
because they violated the essential condition W(0)-0. One way to correct that 
would have been to construct the trial functions more carefully; they could have 
been 


{1 -cos tfx, 1 -cos 2 ttx, 1 -cos 3?rx , . . . } (and also x 2 !) 

The alternative of keeping the extra trial function 1, and imposing W(0)**0 
through a Lagrange multiplier, is equally correct. The example verifies that the 
choice can be based on computational convenience. In other problems it might be 
less easy to adjust the trial functions to satisfy the essential conditions. 

3. The special feature of the norm || W|| 2 - JV^") 2 dx the distance 

from any admissible w to the exact solution W satisfies 

\ II W - w|| 2 = P(w) - P(W) (11) 


It follows that minimizing the potential energy P over all trial functions w 
automatically minimizes the distance to W. If the trial space actually contains 
W, the distance is zero and we have the global minimum of P. In the typical 
case, when the Rayleigh-Ritz method keeps w in a finite -dimensional space, the 
best w is the projection of the exact W onto the trial space. The w that 
minimizes P also has minimum error - but to establish (11) we must measure the 
error |w - w|| in the correct # energy norm'. 

The verification splits off a term of integration by parts: 


w - W|f - 




[(w ,f ) + 2W”(W" - w") - (W") ] dx 
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^ (W - w) dx + [W"(W' 
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w ' ) - W' ' ' (W 


w) ] 
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the boundary terms vanish because of the essential conditions on W and w at x=0 , 
and the natural conditions on W at x-1. Writing p for d*W/dx A , and substituting 
into the first line, we recognize its right hand side as 2P(w) - 2P(W) . The 
energy norm imposes itself. In the < energy Inner product*, there is an 
orthogonal projection of W onto the trial space. That is why Rayleigh-Ritz, and 
finite elements, do so well. 


4. A final small worry. If a set Is not complete, as the original cosines 
were not, the terms they give should satisfy Bessel' s inequality: 


cos nirx|| 2 <||w|| 2 
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This means that the component of W inside the trial space should not be larger 
than W itself. But the inequality was tested in the short section prior to this 
one, and it failed. The left side was -H*> and the integral of (W n ) 2 on the right 
side was finite. After some thought one realizes that the a n may be the cosine 
coefficients of W, but they are not the coefficients in the right inner product! 
Bessel ' s inequality with the series (3) is satisfied in the L 2 norm. In the 
energy norm, we must use the coefficients (7), (8) obtained from the Rayleigh- 
Ritz method. 

Summary 

It follows from the above discussion that the Rayleigh-Ritz method has 
converged upon the best approximation (as measured in the energy norm) in the 
space spanned by the cosines. We can write the approximation (9) in the 
alternate form 


W* 


T, (1 ' cos n * x) 

- n f 


The above is not the cosine expansion of W (eq. 3). The two are different 
because the x 2 term was forgotten. Without that term we have Bessel's inequality 
||w*|| 2 <||w || 2 . With x 2 included to complete the set of trial functions, the 
variational solution will become the exact one: 
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We can verify that p/12 is the Fourier coefficient in the correct inner product: 


(W, x ) JW" • 2dx p 

( x 2, x 2)“ J 2 • 2dx " 2 
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This example serves to illustrate an elusive pitfall in practice. It is possible 
to see rapid numerical convergence and conclude that a good approximation to_ the 
exact solution has been obtained. In reality, if one's set of basis. functions 
are not complete (in the appropriate inner product space) the approximate 
solution may be highly inaccurate. Finite element experts will note that the 
pitfall could have been avoided by applying the patch test. That requires the 
special solution W—x 2 , with constant strain, to be reproduced by the trial 
functions - after imposing conditions on an element boundary consistent with this 
particular W. The cosines would not have reproduced x 2 . Thus the patch test 
would have failed, correctly indicating that the set was Incomplete. 

The reader will have noticed from the beginning that all combinations of the 
cosines satisfy the extra condition W'(1)=0. If that had been the essential 
condition at the right-hand end (see Reference 2, p. 174), with W' ' ' (1)=0 as 
natural condition, the cosines would have been adequate. 
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Note: A completely parallel example can be constructed for the second- order 

equation -u"-l, with essential condition u(0)-0 and natural condition u'(l)=0. 
The corresponding set of trial functions, complete in L 2 but incomplete in Hg 1 , 
is { sin n?rx} . 


Numerical Results 

Here we compare the cosine expansion to a different expansion that the theory 
guarantees to succeed. The latter comes from the eigenfunctions of a uniform 
vibrating free -free beam: 

- /9% with ^"(0) - *'”<0) - ^"(1) - *"’(D - 0 (12) 

dx 

Note that a more natural choice is a clamped-free beam, satisfying * (0)—0 

at the left endpoint (and more like 1 - cos mrx) . We wanted to see how the 
Lagrange multipliers would enforce these geometric conditions, when they are not 
imposed on each eigenfunction. 

The results are striking. For the cosine expansion (through cos Nttx) we 
tabulate the deflections at x=l and their errors (see Table I) . The exact value 
is W*(l)-p/24. W* will be supplemented by the px 2 /12 correction term, which is 
twice as large! Together they reach the correct value W(l)-p/8. 

What is significant is the 1/N 3 convergence rate. Doubling N reduces the 
error in the last column by a factor of 8. 

Contrast that with the results using the free-free eigenfunctions 

<f> x - 1 with - 0 
4> z - * - \ with fi z - 0 

- (sin 0 n - sinh 0 n )(cos 0 n x + cosh /3 n x) 

- (cos p n - cosh £ n )(sin 0 n x + sinh P n x) 
with cos £ n cosh - 1 


The frequencies are asymptotic to 7m + constant. The expansions can almost be 
carried out by hand (with the help of orthogonality), but imposing W(0)«W'(0)«0 
by Lagrange multipliers needs a simple code. The deflections are again tabulated 
at the free end x-=l , and you will notice the change in convergence rate (see 
Table II) . The error decreases like 1/5N. At N— 800 we are far above the error 
achieved previously at N-10. 
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Table I. 


Table II. 


N W* N (l)/p Error N W N (l)/p Error 


10 

0.04165996 

0.00000671 

50 

0.120908 

0.004092 

20 

0.04166582 

0.00000085 

100 

0.122966 

0.002034 

40 

0.04166656 

0.00000011 

200 

0.123988 

0.001012 

80 

0.04166665 

0.00000002 

400 

0.124496 

0.000504 

Exact 

0.04166667 


800 

0.124750 

0.000250 




Exact 

0.125000 



Finite Elements 


This last section looks at the effect of enforcing natural boundary conditions 
(as well as essential conditions) in the finite element method. The trial 
functions will be piecewise polynomials, and for the bending problem the natural 
choice is piecewise cubics. There are two parameters at each node - 
displacement and slope. Because they are continuous between elements, the trial 
functions w are in the class C 1 - with one continuous derivative, and jump 
discontinuities in the second derivative. This guarantees that w will be 
conforming; its bending energy (and therefore its norm!) is finite. We have only 
to think about the boundary conditions. 


The essential conditions fix w 0 -0 and w'g-0 at the left end. With N intervals \ 

of length h-l/N, and two degrees of freedom per node, that leaves a 2N- 
dimensional trial space - provided no conditions are applied at the right end. 

If we do impose the natural boundary conditions, they yield two relations between \ 

w *n-i> w n» w 'n- Therefore the trial space is reduced to dimension 2N >2. 

The question is whether this Is wise - to compel w to meet conditions that we I 

know to be satisfied by the true solution W. { 

A functional analyst would say it is foolish. The smaller trial space (it is : 

a subspace of the larger one) cannot give a better approximation to W. The error I 

must be greater when degrees of freedom are removed. But that is the error In = 

the energy norm, where the identity (11) holds and Rayleigh-Ritz picks the best 
finite element approximation as the projection. We might still hope that 
pointwise, and particularly near x-1 , there is something to be gained by = 

imposing the natural conditions. 

This example is so simple that the calculations can be done with pencil and 
paper. Furthermore, it has the special property of superconvergence . The finite - 

element approximation is exact at the nodes , In the full trial space, with no 
constraints at x-1, there is agreement with the true W at every meshpoint. In 1 

the smaller trial space, the natural boundary conditions require the cubic to be 
a linear function within the final Interval. (Then it satisfies w M =w' f '=0 
throughout the Interval; that Is the price for Imposing those conditions at x=l . ) ~ 

In this case it Is still exact at all other nodes! The discrepancy between the - 

two finite element solutions is small, and very local, but the winner is clear, f 

i 
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Even near x-1, it is better when the natural conditions are left alone. By 
satisfying them, we spoil the accuracy. 

It is a pleasure to verify superconvergence in this case. The element 
stiffness matrix and element force vector are 
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Those refer to the local parameters w ^ , w' w'j. 

the overlapping parts: ^ J 


The assembly combines 


[ i ] r [ [ ] ] 

Therefore the global stiffness equations for j<N (found from 3P/3Wj-=3P/3w' j-0) 
come from the middle rows of that assembly: 


4 

-12wj ^ - 6hw'j ^ + 24Wj - 12Wj +1 + 6hw'j + ^ - ph 
6hWj ^ + 2h^w'j ^ + 8h^w'j - 6hw^ + ^ + 2h^w'j + ^ - 0 


(13) 


Those equations are exactly satisfied by the nodal values Wj and W'j of the true 
solution W in equation (2). 

At the end x-1, imposing w^-0 and w^' '-0 leads to 

W N “ W N-1 and W N “ Vl + hW N-l (14) 


The cubic w is forced to be linear In the last interval. On the other hand, if 
these natural conditions are not imposed, the only difference from (13) is that 
no terms are assembled for the interval beyond the boundary: 


- 12w xt . - 6hw' + 12w - 6hw' - ph^/2 

N-l N- 1 N N 

6hw„ , + 2h 2 w' T , - 6hw + 4h 2 w' ph 5 /12 

N-l N-l N N 


(15) 


The latter are satisfied by the true W but (14) is not - even though W itself 
does satisfy the natural conditions . 

Remark . Even if the right side of W" "-p is not constant in the original 
differential equation, it is still true that the finite element approximation has 
w-W and w'-W' at the nodes. This was observed already by Pin Tong 3 . Those 
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conditions determine one particular interpolating w - and to show that this 
special w is the finite element solution, we verify that its error w - W is 
orthogonal to all trial functions. Then w is the correct projection of W, and 
superconvergence is verified. 

Proof of orthogonality: In each interval the trial functions are cubics C(x), 

and the contribution to the inner product is an integral over that Interval: 


J (w-W) "C"dx - [(w - W)'C" - (w - W)C"'] 
+ J(w - W)C”"dx 


(16) 


this is zero because w-W and (w-W)' are zero at the endpoints, and C"" vanishes 
identically. 

The argument still applies when we impose the natural boundary conditions, and 
force w to be linear in the last interval. Over all other intervals (16) Is 
zero as before. Over the last interval it is zero because C"-0. Therefore the 
cubic which matches W and W' at every node except x-1 is the finite element 
solution in this case. It would have been more accurate to match at every node, 
by not forcing w" (l)-w' ' ' (l)-0. 

In this example, and surely in more realistic and more complicated 
applications of the displacement method, it is better to leave the natural 
boundary conditions to nature. 
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